###################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)

###################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character") ,
    make_option(c("--maf_file"), type = "character") ,
    make_option(c("--geneN"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    input_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.tsv"
    maf_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.maf"
    out_path <-"/public/home/xxf2019/20220915_gastric_multiple/dna_combinePublic/public_ref/images/mutBurden"
    geneN <- "TP53"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
maf_file <- opt$maf_file
geneN <- opt$geneN
mutType <- opt$mutType
out_path <- opt$out_path

dir.create(out_path , recursive = T)

##############################################################################

if(1!=1){
    col <- c(
      brewer.pal(9,"YlGnBu")[6],
      rgb(234,106,79,alpha=255,maxColorValue=255),
      rgb(203,24,30,alpha=255,maxColorValue=255),
      rgb(255,0,0,alpha=255,maxColorValue=255)
      )
}
col <- c(
    brewer.pal(9,"YlGnBu")[6],
    rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
    rgb(red=2,green=100,blue=190,alpha=255,max=255) ,
    rgb(255,0,0,alpha=255,maxColorValue=255)
    )

names(col) <- c("IM" , "Mut" , "Wild" , "GC")

##############################################################################

dat <- data.frame(fread(input_file))
dat_mut <- data.frame(fread(maf_file))

##############################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
dat_mut <- subset( dat_mut , Variant_Classification %in% Variant_Types & Hugo_Symbol == geneN )

##############################################################################

dat_all <- dat
dat_all <- rbind(dat_all , dat)
dat_all <- subset( dat_all , MS_Type == "MSS" )

if(1!=1){
    if(mutType == "Mut"){
        dat_all <- subset( dat_all , Tumor %in% dat_mut$Tumor_Sample_Barcode )
    }else if(mutType == "Wild"){
        dat_all <- subset( dat_all , !(Tumor %in% dat_mut$Tumor_Sample_Barcode) )
    }
}

dat_all$gene_mut <- ifelse( dat_all$Tumor %in% dat_mut$Tumor_Sample_Barcode , "Mut" , "Wild" )

##############################################################################
## 分不同亚型
#dat <- subset( dat , Molecular.subtype %in% c("CIN" , "GS") )
dat <- subset( dat_all , Class!="IGC + DGC" )

##############################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

dat_tmp <- c()
for( classN in unique(unique(dat$Class)) ){

    dat_plot_tmp_use <- subset( dat , Class == classN )

    ###########################################################################################
    ## 样本数量
    #dat_plot_tmp_use$id <- paste0( dat_plot_tmp_use$From , "_" , dat_plot_tmp_use$Class )
    sample_class_num <- dat_plot_tmp_use %>%
    group_by(gene_mut) %>%
    summarize( nums_class = length(unique(Tumor)) )

    dat_plot_tmp_use <- merge( dat_plot_tmp_use , sample_class_num , by = "gene_mut")

    a <- dat_plot_tmp_use[dat_plot_tmp_use$gene_mut==unique(dat_plot_tmp_use$gene_mut)[1],"BurdenExon"]
    b <- dat_plot_tmp_use[dat_plot_tmp_use$gene_mut==unique(dat_plot_tmp_use$gene_mut)[2],"BurdenExon"]
    
    p <- wilcox.test( a , b )$p.value

    if( p < 0.01 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
    }
    
    dat_plot_tmp_use$p_text <- ""
    dat_plot_tmp_use$p_text[1] <- p_text
    dat_plot_tmp_use$Class_use <- paste0( dat_plot_tmp_use$gene_mut , "\n" , "(" , dat_plot_tmp_use$nums_class , ")" )
    dat_tmp <- rbind(dat_tmp , dat_plot_tmp_use)
}

##############################################################################

y_max <- max(dat_tmp$BurdenExon) + 280
dat_tmp$Class_use <- factor( dat_tmp$Class_use , levels = unique(dat_tmp$Class_use)[order(unique(dat_tmp$Class_use) , decreasing=F)] , order = T )
dat_tmp$gene_mut <- factor( dat_tmp$gene_mut , levels = unique(dat_tmp$gene_mut)[order(unique(dat_tmp$gene_mut) , decreasing=F)] , order = T )
y_lab <- "Mutation rate per MB"

plot <- ggplot(dat_tmp, aes(x = Class_use , y = BurdenExon, color = gene_mut , fill = gene_mut)) +
        #geom_boxplot(size = 1.2 , outlier.shape = NA ) + ## 去除散点，加粗线
       # geom_jitter(position = position_jitterdodge(0.8) , size = 1) + 
        geom_violin(trim=FALSE) +
        geom_boxplot(width=0.2,position=position_dodge(0.9),fill="white",color="black")+ #绘制箱线图
        scale_y_continuous(trans='log10', breaks = c(1,10,100), labels = c(1,10,100),limits=c(0.1,y_max)) +
        facet_grid( .~ Class , scales = "free_x" ) +
        scale_color_manual(values=col) +
        scale_fill_manual(values=col) +
        geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
        xlab(NULL) +
        ylab(y_lab)+
        theme_bw() +
        theme(
            legend.position = 'none',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_blank(),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 12,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            axis.text.x = element_text(size = 12,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            strip.text.x = element_text(size = 15, colour = "black",face='bold') ,
            axis.line = element_line(size = 0.5)) 
    
out_name <- paste0(out_path , "/MutationBurden.compare.IGC_DGC.MolecularType.",geneN,".MutWild.pdf")
ggsave(file=out_name,plot=plot,width=7.4/1.5,height=5.4/1.5)